clear

img = 'img5.tif';

if(strcmp(img,'img.tif'))
    c = 1:4;
    t = 1e-3*(1:5:16);
else
    c = 1:5;
    t = 1e-3*(1:5:21);
end

for n = c
    y = imread(img,'Index',n,'Info',imfinfo(img));
    
    % y-signal (in camera frame, not using expt axes)
    yInt = sum(y,2)/1e6;
    xInt = sum(y,1)/1e6;
    
    yC = 1:1040;
    xC = 1:1392;
    
    yMax = find(yInt == max(yInt),1);
    xMax = find(xInt == max(xInt),1);
        
    f = @(x,xdata) x(1)*exp(-(xdata - x(2)).^2 ./ x(3)) + xdata*x(4) + x(5);
    
    ampGuessY = max(yInt) - min(yInt);
    ampGuessX = max(xInt) - min(xInt);
    
    pY = lsqcurvefit(f,[ampGuessY yMax 100 1 min(yInt)],yC',yInt);
    pX = lsqcurvefit(f,[ampGuessX xMax 100 1 min(xInt)],xC,xInt);
    
    wY(n) = pY(3)
    wX(n) = pX(3)
    
    cnX(n) = pX(2);
    cnY(n) = pY(2);
    
    plot(xC,xInt,yC,yInt,xC,f(pX,xC),yC,f(pY,yC))
    pause
end


tt = t.^2;

px2m = (12.9e-6)^2; % 1 pixel = 6.45 um * 2 = 12.90 um. 

wYm2 = wY*px2m;
wXm2 = wX*px2m;

coeffsY = polyfit(tt,wYm2,1);
coeffsX = polyfit(tt,wXm2,1);

v02y = coeffsY(1);
v02x = coeffsX(1);

m = 87*1.66e-27; % mass of Rb87
k = 1.38e-23;

Tx = v02x*m/(2*k)
Ty = v02y*m/(2*k)

r0y = sqrt(coeffsY(2))
r0x = sqrt(coeffsX(2))

plot(tt,wYm2,'.-',tt,wXm2,'.-',tt,polyval(coeffsX,tt),tt,polyval(coeffsY,tt))





